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ABSTRACT  , 

A nev  scalar  labelling  algorithm  is  presented  for  solving  a system 
of  equations  by  simplicial  approximation.  The  method  presented  exhibits 
strong  convergence  behavior  and  supercedes  previous  simplicial  pivot 
algorithms  due  to  the  elimination  of  an  extra  dimension,  the  simplification 
of  the  pivoting  process  by  using  scalar  rather  than  vector  labels,  and,  most 
importantly,  the  nature  of  the  homotopy  path  taken  which  has  the  remarkable 
properties  of  monotonicity  and  Jacobian  invariance.  Examples  are  presented 
vherein  the  nev  method  converges  but  Newton's  method,  Euler's  method  and 
previously  proposed  simplicial  pivot  algorithms  fail  to  converge. 
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1.  INTRODUCTION 


la  this  paper  ve  consider  the  method  of  simplicial  approximation  for 
solving  the  general  problem 

find  x€  satisfying  f(x)  ■ 0 
vhere  f:  Rn  •*  Rn  is  continuously  differentiable. 

A number  of  results  on  simplicial  approximations  appear  in  the  litera- 
ture: [1],  [2],  [U],  [51,  [61,  [7l,  [81,  [91.  [101,  [111,  [12),  [131,  [1U], 

[17],  [18],  [19],  [20].  For  example,  it  is  known  that  a vector  labelling 
[l),  [10]  method  due  to  Merrill  [lU]  tracks  the  "homotopy  path”  f(x)  ■ Xx, 

X <_  0 assuming  (without  loss  of  generality)  that  the  method  is  initiated  at 
the  origin.  In  [9 1 we  showed  that  a "scalar  labelling"  can  be  defined  on  an 
appropriate  triangulation  so  as  to  follow  precisely  the  same  homotopy  path. 

This  is  of  computational  significance  since  it  eliminates  the  need  for  an 
extra  "sandwich"  dimension  and  the  need  to  pivot  on  a linear  system. 

In  [l]  and  [10)  a different  vector  labelling  was  defined  in  such  a way 
that  under  appropriate  assumptions  the  homotopy  path  followed  is  of  the  form 
f(x)  * Xf(0),  0 < X < 1,  assuming  once  again  that  the  starting  point  is  at 
the  origin.  In  this  latter  case  an  efficient  implementation  of  the  algorithm 
suffers  from  a difficulty  in  determining  the  initial  simplex.  However,  once  a 
satisfactory  start  is  obtained,  the  ensuing  path  will  in  general  be  distinctly 
different  than  the  one  generated  by  Merrill's  algorithm  and  it  exhibits  powerful  con- 
vergence behavior.  In  particular,  if  f has  a nonsingular  Jacobian  at  the  origin. 
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and  if  the  origin  is  the  only  preimage  of  f(0),  and,  finally,  if  the  pre- 
image  of  the  line  segment  [0,  f(0)]  (i.e.  the  segment  Joining  0 and  f(0)) 
is  compact  then  the  method  is  assured  to  converge  to  a zero  of  f.  Further- 
more, if  f is  also  1-1  on  the  preimage  of  [0,  f ( 0 ) ] then  the  method 
is  "norm-reducing”  in  the  sense  that  ||f(x)||  decreases  as  one  moves  along 
the  path.  This  is  the  first  known  example  of  a complementary  pivoting 
algorithm  which  exhibits  this  more  classical  monotonicity  behavior.  Another 
remarkable  feature  of  the  algorithm  is  that  the  same  path  is  generated 
regardless  of  whether  the  algorithm  is  applied  to  f(x)  or  to  the  modified 
function  g(x)  obtained  by  premultiplying  f by  the  inverse  of  its  Jacobian 
at  the  starting  point  i.e.  g(x)  » J~^(0)*f(x).  That  is,  the  path  defined 
by  f(x)  » Af(o)  is  identical  to  the  path  defined  by  g(x)  » Ag(0).  This 
invariance  appears  to  be  of  importance,  for  with  most  other  complementary 
pivoting  algorithms  a different  path  is  obtained  by  changing  from  f to 
g,  and  in  fact  the  choice  of  g is  indicated  (see  [8]  and  [9])  because 
of  much  stronger  convergence  behavior.  An  obvious  way  of  stating  essentially 
the  same  point  is  as  follows . Changing  the  subscripts  of  the  f^  functions 
can  change  the  course  of  the  path  f(x)  * Ax,  and  hence  the  subscripting 
can  affect  the  convergence.  The  path  f(x)  ■ Af(0)  is  Independent  of  sub- 
scripting. v 

In  this  paper  we  present  a "scalar  labelling"  approach  that  generates 
the  homotopy  path  f(x)  » Af(o) , 0 <_  A ^ 1,  and  overcomes  the  above  mentioned 

difficulty  in  starting.  The  resulting  algorithm  is  successfully  applied  to 
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two  examples  where  Merrill's  algorithm  appears  to  fail.  This  points  out 
the  advantage  of  following  the  path  of  the  fora  f(x)  * Xf(0),  0 <_  X <_  1. 

Use  of  the  scalar  labelling  to  define  this  path  leads  to  Improved  speed  of 
convergence  because  of  the  ease  in  pivoting  and  the  elimination  of  an 
extra  dimension.  We  also  show  an  example  where  our  algorithm  converges  but 
where  both  Newton' s method  and  Euler's  method  fail. 

2.  DESCRIPTION  OF  THE  LABELLING  AND  THE  INITIAL  SIMPLEX 

Let  f:  Rn  -*•  Rn  be  continuously  differentiable.  Furthermore,  assume 

that  the  following  conditions  hold: 

1)  J-(0)  is  nonsingular  where  J.(0)  is  the  Jacobian 

r r 

of  f at  x ■ 0. 

2)  tJ"1(0)  f(0)]1  #0  all  i - 1,  2 n.  / 

Define  g:  Rn  «*>  Rn  by  the  mapping  I 

g(x)  * J^(0)  f(x)  . / 

Then  by  2)  above  g^(0)  i 0,  all  i.  This  function  g was  first  used  in 
the  complementarity  context  in  [8].  Given  any  fixed  e such  that 
0 < c < min  k|g.(0)|,  define  the  label  of  x to  be 

J * 


[min  {i 


Le(xJ 


g4(x)  6,(x)  .,  &.(*)  ' 

ijTor-5jrirVj}  lf  s^or 6 (0* 1 - forsone 


n + 1 otherwise 


u 

The  change  from  f to  g Is  made  In  order  to  prove  the  existence  of  a 
unique  (n  + 1)  - complete  Initial  simplex  in  a neighborhood  of  the  origin. 

The  following  series  of  Lemmas  will  establish  this  result. 

Let  us  employ  the  definition  | |x| | ■ max  |x. | throughout  the  paper. 

1 

Also,  for  5 > 0,  define  the  6-cube 

D(«S)  « {x  C Rn|-«S  < xt  < 6 , all  1} 

Row,  given  £ > 0,  define 

<jg(e)  * {v®,  v1,  ....  vn}  where  v°  ■ 0,  v1  ■ - e sgn  g^ ( 0 ) e1,  each  1 > 1, 

where  e*  Is  the  i**1  unit  vector  in  Rn. 

We  shall  Impose  a triangulation  on  Rn  In  such  a way  that  aQ(c)  la 
an  n-simplex  of  the  triangulation  and  such  that  each  n-simplex  intersects 
only  a single  orthant  of  Rn.  This  triangulation  is  specified  in  the  appendix. 
Define  a simplex  of  the  triangulation  to  be  (n  + Incomplete  If  the  labels  on 

the  vertices  of  the  simplex  are  1,  2,  . ..  a + 1,  and  define  It  to  be 

n-complete  if  the  labels  are  1,  2,  ....  n.  (Thus,  only  a n-simplex  can  be 
(n  + incomplete,  and  only  (n  - 1)  and  n-slmpllces  can  be  n-complete.) 

We  now  show  that  there  is  a D(6)  such  that  for  all  trlangulatlons 
with  grid  size  small  enough  (as  measured  by  e ) Oq(e)  Is  the  unique 
(n  + l)-complete  simplex  of  D(6).  Recall  that  e is  used  to  define  the 
grid  size  In  the  sense  that  the  vertices  of  the  initial  simplex  cr^fe)  are 
given  by 


/ 
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0 

▼ ■ 0, 


Consider  an  x £ D(6)  for  come  S > 0.  Recalling  that,  by  assumption, 
g^O)  + 0,  all  i,  the  Taylor  expansion  for  gj(x)  about  the  origin  yields 


g.  (x)  x. 

^roT"1*^Toy'R(||x||,• 


i * 1,  2 , « « • , n 


where 


lim 

IxlhO 


Rillxl 


We  choose  6 > 0 in  such  a way  that 


i)  <S  < U/5  |g,(0)| 

i 1 

ii)  for  each  x£D(<S),  ||x||  i 0,  we  have 


Thus  for  each  i. 


R . x I 


Mi^oTT * Mxtr  * “Jn  ^TiJToTT 


-VfcJoTr  * 


Lemma  1 Let  | |x|  | » e £ (0,  <S) , xi  ■ 0.  Then  Lg(x)  + i . 


Proof; 


||^-i*e(iixii) > i - ■ i - “>4 


hence 


Lc(x)  * i 
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Leana  2 Let 
Then 

Proof: 


Lemma  3 Let 
Then 

Proof: 


x € D(S),  IxJ  « | |x|  | 1*0,  and  sgn  xt  « sgn  gt(0). 
L£(x)  i*i  if  0 < e < <5 


8,  (x) 

s[ToT  “ 1 * 


R( 11*11) 


» l ♦ 


I g±  ( o ) I 


♦ H(||x||) 


• 1 + 


1 + 


> 1 » l 1 s{foh  > 1 *»  Le(x)  * 1 if  0<e<5 


x€  D(5),  |x± | » | |x|  | i*  0,  sgn  ■ -sgn  gA(0) . 
L£(x)  i*  n + 1 if  0 < e <,  | |x|  | . 


gJoT 


^■1-t£tW*r(|1 1x111 


- 1 
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...ljML.i.m.i. 


rfipofr 1 • Hl^oTT  - 1 • iTi^oTT 


if  0 < e < xl 


Purtb'r=or'  l^loT  * 1 - t T^oIt*  (4ffloTF  r(||,||)] 

- 1 ■ £ T^ToIr  - 1 • £ T^ToTT  (*ln" 


1*11  < <S) 


> 0 (since  5 < {r  min  | g.  (0 ) | ) - 
5 i 


Thus  L£(x)  i n + 1 


Learn*  o (e)  is  (n  ♦ l)-complete  if  0 < c < <5. 


Proof:  Take  any  i 6 {l,  2,  ....  n}. 

Note  that  |r*|  ■ Hr^H  ■ e and  sgn  ▼*  ■ -sgn  g^(0) . 

Hence,  by  Lemma  3,  L^r1)  f n + 1 

and,  by  Lemme  1,  L£(vi)  f J all  J i i. 

Hence  L^v1)  ■ i.  Since  L(0)  ■ n ♦ 1,  oQ(e)  is  (n  ♦ l)-complete. 


Lemma  5 Given  t £ (0,  <S) , let  O * {u®,  u\  . . . , un}  # O^Ce) , O c D(<S) , 

and  assume  a is  in  the  same  orthant  as  Oq(e).  Then  n + 1^  L£(a), 


Proof:  Take  any  J 6 {0,  1,  ....  n}.  Note  that  ||u^||  + 0. 

Hence  choose  i such  that  |u^|  « ||u^||  f 0. 

Since  a c orthant  aQ(e),  sgn  ■ -sgn  g^fO) 

Since  ||uJ||  >_t,  by  Lemma  3,  LG(u^)  j n + 1 
Thus  n + 1 ^ LG(a)  • 

Next,  consider  an  n-simplex  a ■ {u^,  u1,  un},  oC D(<S)  and'  a 

not  in  the  same  orthant  as  aQ(e),  given  0 < e < 5.  Let  u » (i^,  Ug,  ...»  uq 
he  a vector  defined  hy 


min  u' 


J 


if 


> 0 all  J ■ 0,  1,  .... 


if  uj  < 0 all  J * 0,  1,  ....  n 


Note  that  for  each  i , either  ujJ  0 all  J or  u^  <_  0 all  J . 


Also 


•note  that  for  any  2 vertices  x,  y of  a,  lxjl  ” ly^ I € (-e,  0,  e}  for  all  J. 

Furthermore , for  any  i and  j 

|u^|  is  either  | u.^ | or  | u^ | + e. 


First,  let  us  consider  the  case  where  ||u||  * 0. 
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Lemma  6 Given  0 < e < 6,  let  a * {u°,  u\  ....  un}  d D(  6 ) , 0 not  in 
the  same  orthant  as  Cq(e),  and  | |u| | * 0.  Then  there  is  an 
i £ {l,  2,  ....  n}  such  that  i^’L£(a). 

Proof:  If  | |u| | ■ 0 then  | |u^ | | <_  e for  all  J . 

Since  a orthant  of  0q(e),  there  is  an  i such  that  u^  # 0 implies 
sgn  u^  * sgn  g^O),  for  all  J.  Consider  any  J ■ 0,  1,  ....  n.  If 
| |uJ  1 1 * 0 then  u^  * 0 so  that  L£(0)  * n + 1 j*  i.  Otherwise,  ||u^||  * e. 

If  u^  * 0,  Lemma  1 implies  that  L£(u^)  j*  i.  If  |u^|  * ||u^||  » e,  then 

Lemma  2 implies  L£(u^)  i i*  £21 

Next,  consider  the  case  where  | |u| | + 0. 

Lemma  7 Let  i he  such  that  Ju^  ■ ||u||  *k£  for  some  positive  integer  k. 

Then  for  all  J » 0,  1,  ....  n either  |u||  * ||u^||  or 

lUil  * k ♦ 1 I lu^l  I * (°f  course  ||uJ||  * 0,  since  I luJl  I i I lul  I * 0-) 

Proof:  Take  any  J ■ 0,  1,  ...»  n. 

Then  |u^|  i I |uJ|  | ■ |u£|  (some  k)  < |ujJ  + e <J  |u|  | + c 

Furthermore  |u^|  >_  | u^ | » | |u| | . 

Hence  if  |u^|  < ||uJ||  (i.e.  |u^|  - ||uJ||  - e)  we  have 

|uj|  » I |u| | , | | uJ | | - | |u||  + E 

Thus,  |u^|  ■ kc  and  |ju^||  * (k  + 1)e 

->  |uj|  - l|uJ|| 
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Ltm  8 


Proa 


Lemma  9 


Consider  o * {u^,  u^",  . un},  uc  D(6),  a not  In  the  same 
orthant  as  <Jg(e).  Let  1 he  such  that  |u^|  * | |u|  j j*  0. 
Then  If  sgn  ^ « sgn  g1(0)  ve  have  L£(uJ)  / i all  J * 0 


Take  any  J£(0,  1,  . . . , n} . Since  / 0,  ve  have 
u^  # 0 and  sgn  u^  = sgn  g^(0).  By  Lenina  7 either 

|u^|  * ||uJ||  or  |u^|  » — * j ||uJ||  vhere  k is 
positive  integer  such  that  | |u| | * ke. 


Ca9«  *•  |uj|  - ||uJ|| 

Then  hy  Lemma  -2  L£  ( u*^ ) / 1 . 

csas-fr-  l^|  « lluJ|| 

g,(uJ)  |u^| 

Then  i^or*1*TiiW*!'(|l“  M) 


» 1 ♦ 


k 

k + 1 


R ( | | u^ | | ) (since 


— > i-) 

k + 1 - 2' 


^ 1 since  u^  £ D(3) 
->  le(uj)  i i 


Consider  O as  above.  Let  i be  such  that 

|ujJ  • ||u||  / 0.  Then  if  sgn  u^  * -sgn  g^(0)  ve  have 

L£(u^ ) f n ♦ 1 all  J ■ 0,  1,  . ...  n. 
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Proof: 


Take  any  j€{0,  1,  n}.  Of  course  sgn  u^  » -sgn  g^O) 

By  Lenna  T,  |uJ|*||uJ||  or  luJ I * I |uJ 1 1 for 


some  positive  integer  k. 

J 


Case  a.  | u" | * | | uJ 

Since  ||u^||  ^ e,  Lemma  3 implies  Le(u^)  t n + 1. 

Ca8e  luil  " F+T  I luJl  I * 


gl(uJ) 

^ gjToT  - 1 - 


Since 


Furthermore 


gi(uJ) 

- Tor  -1 


e 


-H-H  I^Tofr  * *B(I|"J|I) 


(since  I |u^ I I < 5) 
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M 


t 


i 1 - k T^foTT  slnoe  SrH  - k 

> 0 (since  iS  < jr  min  |g.(0)|) 

5 i 


->  L£(uj ) + n + 1, 


Theorem  1. 


crg(e)  is  the  unique  (n  + l)-complete  simplex  in 
D(<S)  if  0 < e < <S. 


Proof:  This  follows  from  Lemmas  1*,  5,  6,  8 and  9- 


The  choice  of  e in  Theorem  1 is  crucial  in  the  speed  of  convergence 

of  our  method.  If  e is  too  large,  it  is  possible  that  the  method  cannot  be 

initiated.  If  e is  too  small,  the  method  may  require  a prohibitively  large 

number  of  pivots  to  reach  a terminal  simplex.  Our  test  examples  show  that 

e » min  |g, (0)|  worked  well, 
i 1 


3.  CONVERGENCE 


The  above  Theorem  1 is  all  that  is  needed  to  validate  our  method. 

The  method  is  now  a familiar  one  in  complementarity  theory:  starting  from 

the  simplex  aQ(e)  generate  the  sequence  of  distinct  slmplices  <T1,  a^t  .... 

such  that  0 <y  is  n-complete  for  i ■ 0,  1,  2,  ....  . Terminate 

upon  reaching  for  the  first  time  an  n-simplex  a T which  is  (n  + Incomplete. 

We  may  now  use  any  vertex  of  Om  as  the  next  point  to  restart  the  method. 

▲ 


d. 


r 
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A condition  under  vhich  this  method  is  assured  to  converge  to 
a zero  of  f is  given  by  the  following  Theorem.  First  define  the  set  H 
by 

H ■ {x|f(x)  ■ Xf(0),  0 £ X £ l}  ■ (x|g(x)  ■ Xg(0),  0 £ X £ 1} 

We  define  the  p-neighborhood  of  H to  be  {y j | |y  - x|  | £ p for  some  x £ H}. 

Theorem  2.  Let  f:  Rn  ■»  Rn,  f continuously  differentiable  on  Rn. 

Suppose  f(x)  * f(0)  iff  x ■ 0,  J'^O)  exists, 

[j'^O)  f(0)]i  j*  0 all  i,  and  H is  compact.  Then  there  exists  a solution 
x to  f(x)  * 0 and  furthermore  as  the  mesh  size  goes  to  zero  the  method  will 
converge  to  a solution. 

Proof:  Let  6 > 0 be  such  that  Theorem  1 holds  and  let  be  a 

sequence,  0 < £ 5,  with  decreasing  to  0.  Let  Cj  denote  the  path 

of  slmplices  generated  by  the  method  for  a triangulation  of  mesh  size  e^.  Then 
for  every  p-neighborhood  of  H,  say  R^,  there  exists  a i such  that  i > i 
implies  c Np.  To  see  this  assume  to  the  contrary  that  there  is  a neighbor- 
hood Rp  for  which  fl  (Rn  - Rp)  j*  0 for  infinitely  many  i.  Let 
x*  € 3Rp  n C^.  Then  it  follows  from  the  continuity  of  f on  Rp,  the  nature 
of  the  labeling,  and  the  fact  that  each  simplex  in  each  path  is  n-complete, 
that  every  cluster  point  of  {x*}  is  in  H,  which  is  a contradiction.  It 
follows  from  this  result  that  for  i sufficiently  large  each  path  terminates 
with  a final  simplex,  say  0^ . Since,  by  Theorem  1,  aQ(e)  is  the  unique 


lU 


(n  + l)-complete  simplex  in  the  <5-cube,  each  terminal  simplex  0.  is 


outside  this  cube.  Since  this  terminal  simplex  contains  labels  1 through 

for  all  i or 


6t(x) 


n ♦ 1,  any  x in  the  simplex  must  satisfy  - ^ * 1 

S^x)  _ i 

— ■ 0 for  all  i.  Since  the  unique  preimage  of  g(0)  is,  by  assumption, 

the  origin,  g(x)  is  approximately  zero  for  x in  0^.  Take  any  sequence 
{s*}  such  that  s1  £ 0^.  It  is  clear  that  any  cluster  point  of  {s*}  is 
a z iro  of  f and  g. 


4.  JACOBIAN  INVARIANCE  AND  MONOTONICITY 

It  is  clear  that  finding  a solution  to  the  system 

(1)  f^x)  » 0,  i » 1,  ....  n 
is  equivalent  to  solving  the  system 

(2)  g^x)  ■ 0,  i ■ 1,  ...,  n 

where  g(x)  • J'^Oj’ffx).  Most  simplicial  pivoting  algorithms  follow 
different  paths  depending  on  whether  system(l)or(2)is  being  solved.  Moreover 
it  is  not  unusual  for  the  simplicial  path  associated  with  (1)  to  be  unbounded 
whereas  the  path  for  (2)  converges  (see  [8],  [9],  [13],  [19] » and  [20]).  In 
other  words,  the  transformation  to  system- (2)  is  known  to  improve  convergence. 
For  the  method  presented  in  this  paper  the  simplicial  path  is  the  same, 
regardless  of  vtoether  (l)  or  (2)  is  attacked.  More  precisely,  suppose  the 
labeling  function  L (x)  is  redefined  exclusively  in  terms  of  f rather 


15 
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than  g (i.e.  replace  g^  with  in  the  definition).  If  this  is 

done,  a unique  initial  n + 1-c omplete  simplex  may  not  exist,  and  hence  it 
may  he  either  impossible  or  difficult  to  implement  the  algorithm.  However, 
if  a unique  first  simplex  is  found,  then  the  path  followed  will  coincide 
with  the  path  generated  by  the  method  of  this  paper.  In  other  words,  the 
paths  for  (l)  and  (2)  are  the  same.  Ve  have  defined  the  labeling  function 
in  terms  of  g (i.e.  we  solve  system  (2))  only  to  constructively  prove 
the  existence  of  Og(e).  This  assures  us  of  a start.  The  fact  that  the 
path  is  invariant  under  the  Jacobian  transformation  endows  this  algorithm 
with  a natural  property  not  shared  by  others. 

Another  property  that  the  method  possesses  is  that  if  f is  1-1 
on  the  hoootopy  path  H » {x|f(x)  » Xf(0),  0 < X < l}  then  ||f(x)|| 
decreases  monotonlcally  to  zero  on  that  path.  Monotonlclty  is  assured 
Inasmuch  as  X will  decrease  from  1 to  0 as  one  traces  the  path  H. 

This  property  is  especially  useful  in  restarting  our  method.  In  other 
simpllclal  pivot  techniques,  there  is  no  assurance  that  the  sequence  of 

V 

approximate  solutions,  x (in  the  terminal  simplex  of  the  iteration 
corresponding  to  e^)  will  become  "better”  approximations  to  a true  solution' 
as  k increases.  In  our  method  however,  for  suitably  small,  at  any 
point  x*  in  a terminal  simplex  it  will  be  true  that  1 1 f (x^ ) 1 1 is  less 
than  | |f(x?c'*^)|  | l.e.  in  each  iteration  the  norm  of  f at  a point  in  the 
terminal  simplex  is  less  than  the  norm  of  f at  a point  in  the  initial 
simplex  for  that  iteration.  Thus  each  point  x^  is  an  improved  estimate. 
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5.  EXAMPLES 

Here  we  present  three  examples.  In  each  example  our  method  converged. 
We  also  attempted  to  solve  the  first  two  examples  by  using  Merrill's 
algorithm  [lU]  initiated  at  the  same  starting  point.  In  each  of  these  tvo 
examples  the  variables  on  the  homotopy  path  appeared  to  grow  without  bound 
and  hence  Merrill's  algorithm  terminated  without  convergence.  In  the  third 
example  both  Newton's  method  and  Euler's  method  failed  whereas  ,our  method 
converged. 

Example  I 

f. (x)  ■ xp  - Ex.-  1000  1 < i < n 

1 1 3*1  J ~ “ 

We  used  our  method  on  this  problem  for  n » 10.  The  starting  point 
was  chosen  to  be  x^  ■ 20  all  i.  The  mesh  size  was  chosen  by 

S+l  ’ lgi*xk)l 

vhere  x is  the  approximation  furnished  by  the  previous  iteration.  This 
selection  of  was  the  most  ideal  among  the  options  we  considered.  As 
seen  in  the  examples,  such  a choice  yielded  e^+1  ■ ~ or  better  in  every 
case.  The  results  are  shown  in  Table  I. 

Example  II 

This  example  is  a standard  R(h^)  discretization  of  a two  point 
boundary  problem  (see  Mortf’  and  Cosnard  [15 ] ) 

u"(t)  ■ j (u(t ) + t ♦ l)3  0 < t < 1,  u(0)  ■ u(l)  » 0 . 
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The  resulting  problem  In  the  unknowns 

Vx)  ■ - Vi  - Vi  * 


Is  defined  by 
+ l)3,  1 < k < n 


where  n Is  taken  to  be  10,  and 


n+1 


- 0,  tk  - kh  and  h - ;pry 


The  method  is  started  at  the  point 

x°  • (e, , . .. , e ) where  E.  * t. (t.  - 1) 
i n ill 

and  the  results  are  presented  in  Table  II.  The  system  of  equations  has 

a unique  solution  x*  ■ (e*,  ....  £*)  where  -0.5  < e*  < 0. 

l n — l — 


Example  III 

Let  f,(x)  ■ nx.  + sin  x,  - Ex,  - 100  - i , 1 < i '<  n. 

1 1 2 1 J*i  J ~ ~ 

A solution  (up  to  an  accuracy  of  8.106  * 10”^)  to  this  problem  for  n * 5 1* 
x » (101.043,  101.166,  101.293,  101.423,  101.557).  We  tested  Newton's  method 

xk+1  ■ xk  - J*1^)  f(xk) 

twice  on  the  above  problem  using  initial  points  x°  ■ (20,  20,  20,  20,  20) 
and  x°  - (120,  120,  120,  120,  120). 

In  both  cases,  the  method  failed  to  converge  due  to  some  component  of 


x growing  large . 
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We  also  tested  Euler's  method  of  the  fora 


■ xk  - j;1^)  [f(xki  * !L=-T  «*>)', 


Jc*l  m Jt  _ f(xk) 

where  T is  a given  positive  Integer. 


k ■ 0,  1,  ...»  T - 1 
k * T,  T ♦ 1,  ... 


This  method  is  a "continuation  method"  which  may  be  visualized  as  that  of 
approximating  the  path  f(x)  ■ Af(0),  0 <_  \ £ 1 by  way  of  Newton  directions 

([16],  p.  232).  The  method  converged  from  starting  point  xj  * 120  all  i 
using  T » 100.  However,  it  failed  to  converge  from  x°  * 20  using 
T ■ 100,  1000,  and  10000. 


With  the  same  initial  points,  we  tested  our  method  on  the  problem 
above  first  using  grid  sizes  computed  via 


Vi  * i«i(*k)i 


where  x is  the  approximating  solution  after  malor  iteration  k.  For 
this  approach  our  method  moved  to  points  within  10  units  of  the  solution, 

1r 

but  failed  to  come  closer  because  points  x were  encountered  for  which 


the  preimage  of  g(xK)  is  not  unique. 
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We  then  used  grid  sizes  of  the  fora 

eQ  ■ min  |g. (x°>| 
i 1 

■ ain  min  | «±  ( xk ) | > 

For  both  initial  points  ■ 20  all  i,  « 120  all  i the  method 
converged  to  a solution.  The  results  of  our  tests  are  shovn  in  Tables  III  and 
IV. 


Conclusion 

In  this  paper  ve  presented  a "scalar  labelling"  algorithm  for  solving  a 
system  of  equations  by  way  of  slmpllcial  approximation.  The  method  presented 
possesses  several  remarkable  properties  such  as  monotonicity  and  Jacobian  invari- 
ance. The  algorithm  compares  favorably  with  previous  algorithms  of  Eaves  and 
Saigal  [U]  and  Merrill  [lU]  due  to  the  elimination  of  an  extra  dimension,  the 
simplification  of  thj  pivoting  process  by  using  scalar  rather  than  vector  labels 
and.  most  importantly,  the  nature  of  the  homotopy  path  taken. 


n - 10  x°  « (20,  20,  20,  20,  20,  20,  20,  20,  20,  20) 


ITER. 

SO. 

k 

X 

MESH  SIZE 

. 

l|f<xlt)|| 

NO.  OF  PIVOTS 
TO  REACH  xk 

1 

(8.5,  11. U 11.1*) 

2.8631U 

1*75.369 

758 

2 

(10.68,  10. 3U,  ....  10. 3U) 

1.07039 

296.865 

300 

3 

(10.13,  10. 3U,  ....  10. 3U) 

.187178 

5»*.  3887 

20 

k 

(10.30,  10. 3U 10. 3U) 

.088 

15.1*807 

10 

5 

(10.28,  10.30,  ....  10.30) 

.020 

6.3 

1*58 

6 

(10.29,  10.30,  ....  10.30) 

.0085 

2.67 

10 

T 

(10.29,  10.3,  ....  10.3) 

7.95  * 10'5 

.021*6 

210 

Table  1 


I|f(xk)|| 


4.56  x IO-2  .1019 

2.24  x io"2  3.26  x io‘ 
1.115  * 10“2  2.25  x 10' 

5-56  x io“3  8. 

2.77  x IO-3  3. 

1.38  x io'3  3. 

6.926  x i(fU  i. 

f2,  -8.2  x io“2,  -.11 
Table 


n * 5 


x°  « (20,  20,  20,  20,  20) 


ITER. 

NO. 

l|f(xk)|| 

*k 

NO.  OF  ITERS. 
TO  REACH  xk 

1 

U1.2U 

1*0 .2U 

95 

2 

6.161 

1.2575 

1028 

3 

6.923 

.629 

1*6 

U 

U.105 

.311*3 

75 

5 

3.691 

.1572 

TT 

6 

3.08U 

.0786 

9U 

7 

l.OUU 

.0393 

115 

8 

.U968 

.0196 

72 

9 

.8970 

9.82  * 10-3 

112 

10 

.1*809 

1* . 91  * 10"3 

213 

11 

.3867 

2.1*6  * 10"3 

271* 

x11  • (101.095,  101.166,  101.291,  101. U21,  101.556) 


Table  III 


x°  « (120,  120,  120,  120,  120) 


ITER. 

BO. 

NO.  OF  ITERS. 
TO  REACH  xk 

1 

2.8ll»62 

6.3276 

95 

2 

25.1*632 

1.5819 

1*3 

3 

5.3331 

.79095 

28 

U 

3.3965 

.3955 

36 

5 

1.3965 

.1977 

51 

6 

.50790 

.0989 

1*1 

7 

.1*579 

.01*91* 

1*8 

8 

.2582 

.021*7 

80 

9 

.0510 

.012U 

68 

x9  - (101. 0W6,  101.165,  101.3,  101. *21,  101.55) 


Table  XY 


The  pivot  rules  are  given  as  follows.  Suppose  {u^,  u1,  u11} 

and  (s, , ....  s ) are  given.  If  the  vertex  u*  is  dropped  we  must 
l n 

specify  a new  initial  vertex  u^  and  a new  permutation  (s^  s2>  ...»  sq). 

The  rule  u*+1  * u*  + will  then  determine  the  new  simplex  (u°,  u\  ...»  un). 

Case  1:  Drop  u1,  l<i  <n-l 

u°  » u°,  st  * s1+1,  si+1  ■ s1,  Sj  ■ Sj  for  J ^ {i,  i + 1}  . 

Case  2:  Drop  u° 

aO  i A « ^ _ A 

U " U * aJ  * 3J+1‘  X1  J in  - 1,  sn  - Si 

Case  3 : Drop  un 

G°  - u°  - Q(s  ),  s * s , s - s , 2 < J < n . 

n 1 n j j-1  — — 

Hote  that  CTg(e)  is  determined  by  u^  * 0 and  the  permutation  (l,  2,  ...»  n) 
i.e.  aQ(e)  * {v°»  v1,  • .-i  v11).  v°  * 0,  y1  ■ - e sgn  g1(0)  e1,  1 <_  i <_  n. 


Appendix:  The  Triangulation  and  the  Pivot  Rules 


The  triangulation  and  pivot  rules  are  given  for  Rn.  Let  Q denote  the 
n x n matrix  whose  i**1  column  is  denoted  by  Q(i),  where  Q(i)  ■ y1  - y*"*\  and 
where  y°  * 0,  v*  ■ - e sgn  g^O)  e*,  ' 1 < i £ n.  That  is 


Let  (s,  , ....  s ) denote  a permutation  of  the  integers  (l,  2,  ....  n).  That 
1 n 

is,  each  st  £ {l,  ...»  n}  and  si  ■ s^  <**>  i * J.  Now  consider  a triangulation 
of  mesh  e.  And  let  {u°,  ....  u11}  be  an  n-simplex  in  the  triangulation, 
where  u°  is  any  point  in  the  e-lattice  (each  is  an  integral  multiple 

of  e).  Associated  with  this  siaiplex  is  a permutation  (s^,  Sg,  ...,  sq)  such 
that 

ui+1  - u1  + <3(s1+1).  0 < i < n - 1. 

That  is,  any  simplex  in  the  triangulation  is  defined  by  the  initial  yertex 

u°  and  the  permutation  (s, , ...,  s ). 

i n 


t 
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